function [reg_lin,mis,reg_low,reg_up,reg_asym,...
    reg_low_size,reg_up_size,reg_asym_size,...
    reg_low_cap,reg_up_cap,reg_asym_cap,...
    reg_low_all,reg_up_all,reg_asym_all,reg_firms,...
    reg_lin1,mis1,reg_low1,reg_up1,reg_asym1,...
    reg_low_size1,reg_up_size1,reg_asym_size1,...
    reg_low_cap1,reg_up_cap1,reg_asym_cap1,...
    reg_low_all1,reg_up_all1,reg_asym_all1,...
    reg_lin2,mis2,reg_low2,reg_up2,reg_asym2,...
    reg_low_size2,reg_up_size2,reg_asym_size2,...
    reg_low_cap2,reg_up_cap2,reg_asym_cap2,...
    reg_low_all2,reg_up_all2,reg_asym_all2]=...
    MAIN_ESTIMATION_EFFICIENT(or,oecd,annual,ed,in,r,quart,...
    ham,hp,fi,siz,ni,ser,cogs,markup,ind,nonlin,near_var,time_fe,...
    chol,bcdz,jq,cmr,mpt,controls,one_step,labor,size_subset,rollover_subset,ww,rep,fe)
gls=0;
% annual=1;%if it's 1 then we ctot model
err_cor=0;%if it's 1 then correlation between two errs is computed
lev_state=0;%use leverage as state;
control_leverage=0;%if 1 put levergae as control but also amy other qarterly variable from data2
joint_two=0;%if joint_two=1 put co12 10 and code will load other gdpp ppp variables
extract_gdp_ppp=0;%if 1 then extracts all required gdp ppp variables for joint_two estimation; MAKE SURE TO PUT CO12=57
if joint_two==1
    extract_gdp_ppp=0;
end

treat=0;%if 1 then include teatment
joint=0;%if 1 include ppp gdp
quinn=0;%use quinn measure
jay=0;
if jay==1
    co1=1;%33 is hard peg, 34 is soft peg, 40 is neither i.e. strongly felxible against either soft or hard peg
else
    if oecd==1
co1=[1];%determines baseline state variable in joint estimation co1=10 takes short rate

    else
    co1=1;
    end
    end
co12=[1];%determines what other state variable should be included in joint estimation;10 is cc;57 is gdp ppp;

if lev_state==1
    co1=1;
end
if extract_gdp_ppp==1;
    co12=57;
end
if joint_two==1;
    co12=10;
end
if control_leverage==1;
    co12=57;
end

no_interaction=0;
if  any(or==[21:23,38])
month=1;  in=1;
   ed=1;%41 end in 2007q4;put 1 for baseline or 41 if you want pre_2008_fundamental results 
else
    month=0;
end
% if quinn==0
%FOR DATA_EU_US: 1 IS NFC SPREAD; 2 IS NFC RATE; 3 IS UNEMP; 4 IS CPI;
% 5 IS SHORT RATE; 6 IS LONG RATE; 7 IS BROAD EFFECTIVE ER; 
%8 IS REAL BROAD EFFECTVE ER; 9 IS EXPORTS; 10 IS IMPORTS; 11 IS SP 
%12 is high yield; 13 is 10-year inflation expectations 
%14 is real 10-yr rate; 15 is nominal 10 yr-rate; 16 is U.S. inflation
%indexed yield; 17 is loans (MFI FOR EA BANKS FOR US);
%18 is like 17 but excludes reverse repos for EA; 19 is IP; 20 is tm gz EA EBP
%21 is asw EA ebp; 22 is cs; 23 is ytm ebp controlling also for macro risk
if oecd==1  
% load data_oecd;%or=94 is offical rr er or=95 is parallel rr data
% % or=2;%1 is cpi, 2 is gdp, 3 is cons, 4 is I, 5 is tb/gdp, 6 is real er, 7 is nominal er
%     r=[9,10,16,32:35];%FOR BASELINE TAKING IT UKRAINE (42) DOESNT KILL RESULTS 
if annual==0
    if ww==0
    if ni==0&&ser==0&&cogs==0&&markup==0&&ind==0&&nonlin==0&&near_var==0&&controls==0&&labor==0&&size_subset==0&&rollover_subset==0
load data_cs_tfp;%or=94 is offical rr er or=95 is parallel rr data
    end
    if ni==1&&ser==0&&cogs==0&&markup==0&&ind==0&&nonlin==0&&near_var==0&&controls==0&&labor==0&&size_subset==0&&rollover_subset==0
load data_cs_tfp_ni;%or=94 is offical rr er or=95 is parallel rr data
    end 
     if ni==0&&ser==1&&cogs==0&&markup==0&&ind==0&&nonlin==0&&near_var==0&&controls==0&&labor==0&&size_subset==0&&rollover_subset==0
load data_cs_tfp_ser;%or=94 is offical rr er or=95 is parallel rr data
     end 
     if ni==0&&ser==0&&cogs==1&&markup==0&&ind==0&&nonlin==0&&near_var==0&&controls==0&&labor==0&&size_subset==0&&rollover_subset==0
load data_cs_tfp_cogs;%or=94 is offical rr er or=95 is parallel rr data
     end 
    if ni==0&&ser==0&&cogs==0&&markup==1&&ind==0&&nonlin==0&&near_var==0&&controls==0&&labor==0&&size_subset==0&&rollover_subset==0
load data_cs_tfp_markup;%or=94 is offical rr er or=95 is parallel rr data
    end 
    if ni==0&&ser==0&&cogs==0&&markup==0&&ind==1&&nonlin==0&&near_var==0&&controls==0&&labor==0&&size_subset==0&&rollover_subset==0
        load data_cs_tfp_ind_def;%or=94 is offical rr er or=95 is parallel rr data
    end 
    if ni==0&&ser==0&&cogs==0&&markup==0&&ind==0&&nonlin==0&&near_var==0&&controls==1&&labor==0&&size_subset==0&&rollover_subset==0
load data_cs_tfp_with_leverage;%or=94 is offical rr er or=95 is parallel rr data
    end
    if ni==0&&ser==0&&cogs==0&&markup==0&&ind==0&&nonlin==0&&near_var==0&&controls==0&&labor==1&&size_subset==0&&rollover_subset==0
load data_cs_tfp_labor;%or=94 is offical rr er or=95 is parallel rr data
    end 
    if ni==0&&ser==0&&cogs==0&&markup==0&&ind==0&&nonlin==0&&near_var==0&&controls==0&&labor==0&&size_subset==1&&rollover_subset==0
load data_cs_tfp_size_subsample;%or=94 is offical rr er or=95 is parallel rr data
    end 
    if ni==0&&ser==0&&cogs==0&&markup==0&&ind==0&&nonlin==0&&near_var==0&&controls==0&&labor==0&&size_subset==0&&rollover_subset==1
load data_cs_tfp_rollover_subsample;%or=94 is offical rr er or=95 is parallel rr data
    end 
    if (ni==0&&ser==0&&cogs==0&&markup==0&&ind==0&&nonlin==0&&near_var==1&&controls==0&&labor==0&&size_subset==0&&rollover_subset==0)||one_step==1
load data_cs_tfp;%or=94 is offical rr er or=95 is parallel rr data
       parfor i=1:size(data2,3)
    if       size(find(~isnan(data2(:,1,i))),1)>80
        tu(i)=i;
    else
                tu(i)=0;
        end
       end
       end
    else
    load data_cs_tfp_ww;
    end
if nonlin==1
      load data_cs_tfp;
end
if labor==1
      load data_cs_tfp_labor;
end  
    
    r=[];%WHAT FIRMS TO OMIT
    if (near_var==1||one_step==1)&&labor==0
%         r=[find(tu==0),43,90,232,384,446,827,941,1034,1123,1195,1255,...
%             121,191,207,697,934,1115,107,178,206,930,1023,1036,1075,1111,1181,1240];
% r=[find(tu==0)];
r=find(tu==0);

   %THE OTHER OMISSIONS ARE FOR FOR FIRMS WHOSE 21ST HORIZON RESPONSE
   %EXCEEDS 0.01 IN ABSOLUTE VALUE
    end
    if labor==0
   if mpt==1
    parfor i=1:size(data2,3)
    if       size(find(~isnan(data2(1:end-19,1,i))),1)>40
        tu(i)=i;
    else
                tu(i)=0;
        end
    end
    r=find(tu==0);
   end
   if cmr==1
    parfor i=1:size(data2,3)
    if       size(find(~isnan(data2(1:end-29,1,i))),1)>40
        tu(i)=i;
    else
                tu(i)=0;
        end
    end
    r=find(tu==0);
   end
   if bcdz==1
    parfor i=1:size(data2,3)
    if       size(find(~isnan(data2(1:end-27,1,i))),1)>40
        tu(i)=i;
    else
                tu(i)=0;
        end
    end
    r=find(tu==0);
   end
   if jq==1
    parfor i=1:size(data2,3)
    if       size(find(~isnan(data2(1:end-29,1,i))),1)>40
        tu(i)=i;
    else
                tu(i)=0;
        end
    end
    r=find(tu==0);
   end
    end
%     if one_step==1&&labor==1
%         r=[1:4];
%     end
mean_mpk_firm(:,r)=[];data2(:,:,r)=[];data3(:,r)=[];
cap_share(:,r)=[];size1(:,r)=[];
[mpk, mpk_ord]=sort(mean_mpk_firm);
[~, size_ord]=sort(size1);
[~, cap_share_ord]=sort(cap_share);
if ww==0
[~,all_ord]=sort(data3);
else
 [~,all_ord]=sort(ww_index);
end
% if labor==0
% data2=data2(:,:,mpk_ord);
% data3=data3(:,mpk_ord);
% end
    omit=0;
    if month==1
        hor=60;
    else
        hor=20;
    end
else
    load data_eme_annual;
%     r=[47];
    omit=0;hor=5;
end
else
    load data_ea_us;
%     r=[];
    omit=0;%omit certain percentile by explicitely defining 1-fu
hor=60;
end
% if time_fe==1
%     for i=1:size(si,2)
%      y=num2str(si(i));
% si(:,i)=str2num(y(1:3));%TRANSFORM INTO 3-DIGIT INDUSTRY (1:3 is three digits, 1:2 if two digit, etc.)
%  end
if labor==0
    si(:,r)=[];
       [B,I] = sort(si);
%        if one_step==0
[~, w] = unique( B, 'stable' ); w=[w;size(B,2)+1];
%        else
%     w=cumsum(ones(size(data3,2)+1,1));
   data2=data2(:,:,I);data3=data3(:,I);
  
%        end
else
%    if one_step==0
       w=[1;7];
%    else
%     w=cumsum(ones(size(data3,2)+1,1));   
%    end
    
end

cc_dep=0;%put cc variable as dependent variable
if controls==0&&one_step==0
nlag=1;%1+NUMBER OF LAGS
else
    nlag=5;
end
% ed=37;
eq=0;%check equity controls
% in=1;%initial starting observation
% hp=0;%take hp
% ham=0;%take hamilton filter
% quart=1;%polynomial detrending
trend=0;%add trend in estimation itself
cont=0;%continuous interaction
% ed=1;
if cont==1
    ramey=0;
end
lev=0;
if oecd==1
    if annual==0
if or~=16
lev=0;%1 implies levels specification
else
    lev=1;
end
    end
else
    if or==1||or==13||or==20||or==21||or==22||or==23
    lev=1;
    else
        lev=0;
    end
end

gz=1;%use gz shock

% if  any(or==[17:18])
%     data2([1:118,1:52,1:98,1:92],or,[1,3,5,15])=nan;
% end
data1=data1(1:end-ed,:,:);
data2=data2(1:end-ed,:,:);
% if or==26||or==28||or==27||or==40
% data2(:,or,[14,25,37])=data2(:,or,[14,25,37])/1000;
% end

data2=data2(1:end,:,:);%EQUATES OUTCOME DATA SIZE WITH SHOCK'S
t=reshape(data2(in:end,or,:),length(data2(in:end,or,1)),size(data2,3));%rows are gdp over time for each cross section given by the columns

%  t1=reshape(data2(in:end,12,:),length(data2(in:end,12,1)),length(data2(1,12,:)));%rows are gdp over time for each cross section given by the columns
if annual==0
if nonlin==0&&chol==0&&bcdz==0&&jq==0&&cmr==0&&mpt==0&&one_step==0
  load data_baseline_with_rpi;
% load data_baseline_with_rpi_solow;
data(:,1)=[data(1,1);diff(data(:,1))];
data(:,3)=[];
% data=data(:,[3,4,5,7,1,2,6,8,9]);
[~,~,usa]=bootstrap_ISTC_news(data(1:end,:),1,5,2,3,1,3,1,1,1,1,size(data,2)-1,0);
t1=usa.shocks(:,1);
t1=repmat([nan(5,1);t1],1,size(data2,3));
  end
    if nonlin==1&&chol==0&&bcdz==0&&jq==0&&cmr==0&&mpt==0&&one_step==0

t1=VAR_SHOCKS_NONLIN(0,1,0,0);
t1=repmat([nan(5,1);t1],1,size(data2,3));  
end
 if nonlin==0&&chol==1&&bcdz==0&&jq==0&&cmr==0&&mpt==0&&one_step==0
  load data_baseline_with_rpi;
data(:,1)=[data(1,1);diff(data(:,1))];
data(:,3)=[];data=data(:,[3,4,5,7,1,2,6,8,9]);
[~,~,usa]=bootstrap_ISTC_news(data(1:end,:),1,5,2,3,1,3,1,1,1,1,size(data,2)-1,0);
t1=usa.shocks(:,5);
t1=repmat([nan(5,1);t1],1,size(data2,3));
  end 
if nonlin==0&&chol==0&&bcdz==1&&jq==0&&cmr==0&&mpt==0&&one_step==0
  load data_baseline_with_rpi_bcdz;
  data(:,1)=zscore([data(1,1);diff(data(:,1))]);data(:,3)=[];
[~,~,usa]=bootstrap_ISTC_news(data(1:end,:),1,5,2,3,1,3,1,1,1,1,size(data,2)-1,0);
t1=usa.shocks(:,1);
t1=repmat([nan(81,1);t1;nan(27,1)],1,size(data2,3));
  end
 if nonlin==0&&chol==0&&bcdz==0&&jq==1&&cmr==0&&mpt==0&&one_step==0
  load data_baseline_with_rpi_jq;
  data(:,1)=zscore([data(1,1);diff(data(:,1))]);data(:,3)=[];
[~,~,usa]=bootstrap_ISTC_news(data(1:end,:),1,5,2,3,1,3,1,1,1,1,size(data,2)-1,0);
t1=-usa.shocks(:,1);
t1=repmat([nan(50,1);t1;nan(29,1)],1,size(data2,3));
  end 
if nonlin==0&&chol==0&&bcdz==0&&jq==0&&cmr==1&&mpt==0&&one_step==0
  load data_baseline_with_rpi_cmr;
  data(:,1)=zscore([data(1,1);diff(data(:,1))]);data(:,3)=[];
[~,~,usa]=bootstrap_ISTC_news(data(1:end,:),1,5,2,3,1,3,1,1,1,1,size(data,2)-1,0);
t1=usa.shocks(:,1);
t1=repmat([nan(37,1);t1;nan(29,1)],1,size(data2,3));
end 
  if nonlin==0&&chol==0&&bcdz==0&&jq==0&&cmr==0&&mpt==1&&one_step==0
  load data_baseline_with_rpi_textual;
  data(:,1)=zscore([data(1,1);diff(data(:,1))]);data(:,3)=[];
[~,~,usa]=bootstrap_ISTC_news(data(1:end,:),1,5,2,3,1,3,1,1,1,1,size(data,2)-1,0);
t1=usa.shocks(:,1);
t1=repmat([nan(33,1);t1;nan(19,1)],1,size(data2,3));
  end
    if nonlin==0&&chol==0&&bcdz==0&&jq==0&&cmr==0&&mpt==0&&one_step==1
 load data_baseline_with_rpi;
t1=[data(1,1);diff(data(:,1))];
t1=repmat([nan(1,1);t1],1,size(data2,3));
   
    end
%   t1=t1(nlag:end,:);
if month==1
    load ebp_sp_monthly;
    %NEXT ROW IS FOR FEBP; ONE AFTER THAT IS IP; ONE AFTER THAN IS
    %UNEMPLOYMENT
    if or==38
     z=ebp;lev=1;quart=0;
    end
    if or==21
     z=ebp_f;lev=1;quart=0;
    end
    if or==22
     z=ip;lev=0;quart=0;z(isnan(ebp_f))=nan;
    end
    if or==23
     z=un;lev=0;quart=0;z(isnan(ebp_f))=nan;
    end
    
z(any(isnan(z), 2), :) = [];
if or==38
 z([end-14:end,1:12],:)=[];
 % IF YOU WANT PRE_2008 RESULTS UNCOMMENT LINE BELOW (AND PUT ED=41
% ABOVE IN LINE 41
%  z(end-119:end,:)=[];%put 120 for pre_2008_fundamental results case
 %and in tandem put ed=120 in main code FIGURE_GENERATE_PRE_2008_FBP_EBP
end 

t=repmat(z,1,42);
end
else
 t1=-VAR_SHOCKS_ANNUAL(oecd);t1=t1(in:end,:);
end

% if or==12
% t=(t-repmat(nanmean(t),size(t,1),1))./repmat(nanstd(t),size(t,1),1);
% t=exp(t)./(1+exp(t));
% end
% t(:,r)=[];data3(:,r)=[];

if annual==0
% t1(:,r)=[];
end
% t1(:,r)=[];%I take 47 out lready in the shock extraction code
% data2(:,:,r)=[];
if hp==1
    for i=1:size(t,2)
          z=t(:,i);
            select = ~isnan(z) ;y = z(~isnan(z));
            if all(isnan(z))
                t(:,i)=nan;
            else
         y=y-hpfilter(y,1600);
         h(i)= adftest(y);

%                   y=y-ONE_SIDED_HP(y,1600);

            end
z(select)=y;
t(:,i)=z;
         
    end
    share_stationary=sum(h,2)./size(data2,3);
end
if ham==1
for i=1:size(t,2)
          z=t(:,i);
            select = ~isnan(z) ;y = z(~isnan(z));
            if all(isnan(z))
                t(:,i)=nan;
            else
                p=4;h=16;
u=z(p+1+h:end,:); 
for j=1:p
uu(:,j,i)=z(p+1-j:end-j-h,:);
% f=uu(:,j,i);
% f(isnan(z(1:end-5,:)))=0;
% f(isnan(z(5:end,:)))=0;
% 
% uu(:,j,i)=f;
end
c=ones(size(u,1),1);
% c(isnan(z(1:end-5,:)))=0;
% c(isnan(z(5:end,:)))=0;

% c=[c cumsum(c)];
x1=[c uu(:,:,i)];
% x1=x(1:end-h,:);
%   x1(isnan(repmat(u,1,size(x1,2))))=0;
% x1(isnan(x1))=0;
% u(isnan(u))=0;u(isnan(x1(:,1)))=0;
% x=[nan(p,size(x,2));x];u=[nan(h+p,1);u];
[~,~,y]=regress(u,x1);
% y(y==0)=nan;
 y=[nan(h+p,1);y];
t(:,i)=y;

            end

%  z(select)=y;
% t(:,i)=y;
%   if all(isnan(z))
%                 t(:,i)=nan;
%   end
 end
end  
 % [r rr]=find(isnan(t(nlag+1:end,:)));
% [r4 rr4]=find(isnan(t(2:end-nlag+1,:)));
%t=t-repmat(nanmean(t,1),length(t(:,1)),1);
if cc_dep==1
    load cc_cubic
    t=cc;
end
% tt=t;
 tt=t;
% load res
if quart>0
% if or<50&&or>40&&or~=47
%     data2(:,or,:)=log(data2(:,or,:));
% end
bb=reshape(data2(in:end,or,:),size(data2(in:end,or,:),1),size(data2,3));
c=ones(size(data2(in:end,or,:),1),size(data2,3));cc=cumsum(c);
if quart==1
parfor i=1:size(data2,3)
       [~,byyy7(:,:,i),res(:,i)]=regress(bb(:,i),[c(:,i) cc(:,i) cc(:,i).^2 cc(:,i).^3]);
end
clc
end
if quart==2
for i=1:size(data2,3)
       [~,byyy7(:,:,i),res(:,i)]=regress(bb(:,i),[c(:,i) cc(:,i) cc(:,i).^2]);
end
end
tt=res;
end
sample=sum(sum(~isnan(tt),2));
tt(isinf(tt))=nan;
t(isinf(t))=nan;
% save('tt','tt');
y=t(nlag+1:end,:);%dependent variable
pp=t(1:end-nlag,:);pp1=t(nlag+1:end,:);
% t1(1,:)=0;
% if lev==0
pp2=t1(2:end-nlag+1,:);
% else
% pp2=t1(1:end-1,:);
% end

pp3=t1(nlag+1:end,:);

 y(isnan(pp))=0;y(isnan(pp1))=0;
 y(isnan(pp2))=0;
 y(isnan(pp3))=0;
 t(isnan(t))=0;
 t1(isnan(t1))=0;
% lev=0;
if one_step==1
       load data_baseline_with_rpi;
% load data_baseline_with_rpi_solow;
data(:,1)=[data(1,1);diff(data(:,1))];
data(:,3)=[];data(:,end)=[];data=[nan(1,size(data,2));data];
data=repmat(data,1,1,size(data2,3));
end
if lev==0
for i=1:nlag-1
    if controls==0&&one_step==0
y1(:,:,i)=diff(t(nlag-i:end-i,:));%lagged dependent variable in first differences
    else
        if controls==1
 y1(:,:,i)=[reshape(data2(nlag-i+1:end-i,6,:),size(data2,1)-nlag,...
     size(data2,3)) reshape(diff(log(data2(nlag-i:end-i,2,:))),size(data2,1)-nlag,size(data2,3))];%lagged dependent variable in first differences
        end
             if one_step==1
 y1(:,:,i)=[reshape(data(nlag-i+1:end-i,:,:),size(data2,1)-nlag,...
     size(data2,3)*size(data,2))];%lagged dependent variable in first differences
        end
  
        end
     end
  else
   for i=1:nlag-1
y1(:,:,i)=t(nlag-i+1:end-i,:);%lagged dependent variable in first differences
   end
end



% d=reshape(data2(1:end,co1,:),length(data2(1:end,1,1)),length(data2(1,1,:)));%ERR data;

   

if month==0
    cc=reshape(data2(in:end,co1,:),length(data2(in:end,1,1)),length(data2(1,1,:)));%capital control data;
cc=ones(size(cc));
else
 cc=ones(size(t));
 end
% cc(:,r)=[];


    cc1=cc;
pp5=cc(nlag+1:end,:);

y(isnan(pp2))=0;y1(isnan(pp2))=0;y1(isnan(y1))=0;

% ccc=cc(1:end-nlag,:);ccc1=cc1(1:end-nlag,:);
% cc_final=reshape(cc,numel(cc),1);cc_final1=reshape(cc1,numel(cc1),1);

% threshold=quantile(cc_final,qu);threshold1=quantile(cc_final1,qu1);
% threshold2=quantile(cc_final1,qu2);threshold3=quantile(cc_final,qu3);
% rrr=find(cc<0.5); 
% 
% 
%    rrr1=find(cc<threshold3);
    
fu=zeros(size(cc));fu11=zeros(size(cc));

% fu11=zeros(size(cc1));
% fu22=zeros(size(cc1));
%  fu(:,[10])=1; % This is the dummy variable that defines the two states
%DO ZLB
  fu(cc<5)=1;%ZLB
% fu(:,14)=1;%JAPAN
% fu(66:end,[1:3,7:9,12:13,16,18,22,23])=1;%ECB
% fu(:,9)=1;%GREECE
%  fu(:,[5,7,9,10,11,13,15,22:26])=1;%fixed EMEs
% fu(:,26)=1;%Turkey
% fu(:,[172:173,83,59,63])=1;%MEXICO
%  fu(:,end)=1;%US

 fu=fu(1:end-nlag,:);fu1=fu;fu111=fu;fu222=fu;
%  fu11=fu;
%   fu11(66:end,[1:3,7:9,12:13,16,18,22,23])=1;%ECB for omit=1
if oecd==1
 fu11(cc>4.9)=1;%floating EMEs
fu11=fu11(1:end-nlag,:);
 fu22=fu;
else
    fu11(:,1)=1;fu11=fu11(1:end-nlag,:);
 fu22=fu;
end

fu(isnan(pp))=0;fu(isnan(pp2))=0;fu(isnan(pp3))=0;
fu(isnan(pp1))=0;fu11(isnan(pp))=0;fu11(isnan(pp2))=0;
fu11(isnan(pp1))=0;fu22(isnan(pp))=0;fu22(isnan(pp2))=0;
fu22(isnan(pp1))=0;fu11(isnan(fu11))=0;fu(isnan(fu))=0;
if joint==1
    fu(isnan(pp3))=0;fu11(isnan(pp3))=0;fu22(isnan(pp3))=0;
end
tot_obs_cc=sum(sum(fu));
tot_obs_cc_share=tot_obs_cc/nnz(y);
tot_obs_cc1=sum(sum(fu11));
tot_obs_cc_share1=tot_obs_cc1/nnz(y);
total_sample=nnz(y);
total_sample_extended=nnz(t);

%fu(r4,rr4)=0;
% tot_obs_cc=sum(sum(fu));%fu is the constant;

for i=1:nlag-1
    fu1(:,:,i)=fu;%fu1 multiplier lagged outcome variable
        fu111(:,:,i)=fu11;%fu1 multiplier lagged outcome variable
        fu222(:,:,i)=fu22;%fu1 multiplier lagged outcome variable
        if joint==1||eq==1||gls==1
               fu1_14_15(:,:,i)=fu_14_15;%fu1 multiplier lagged outcome variable
         end
end
for i=1:nlag
    fu2(:,:,i)=fu;%fu2 multiplies gz
        fu112(:,:,i)=fu11;%fu2 multiplies gz
                fu1122(:,:,i)=fu22;%fu2 multiplies gz
                if joint==1||eq==1||gls==1
               fu2_14_15(:,:,i)=fu_14_15;%fu1 multiplier lagged outcome variable
         end
end
a=ones(size(y,1),size(y,2));
a(isnan(pp))=0;a(isnan(pp2))=0;a(isnan(pp1))=0;a(isnan(pp3))=0;
  if joint==1
      a(isnan(pp3))=0;
  end
  a=reshape(a,size(y,1),1,size(y,2));

a1=ones(size(y,1),size(y,2));a1(isnan(pp))=0;
a1(isnan(pp1))=0;a1(isnan(pp2))=0;a1(isnan(pp3))=0;
c2=fu.*a1;
if omit==0
c1=(1-fu).*a1;
else
   c1=fu11.*a1;
end
 if trend==1
 c22=fu11.*cumsum(a1);%TREND
 else
c22=fu11.*(a1);
 end
% if time_fe==1
%     u=zeros(size(y,1),1);
%   for i=1:size(y,1)
%       c22(:,:,i)=repmat(u,1,size(y,2));
%       c22(i,:,i)=1;
%   end
% else
%  c22=fu11.*(a1);   
% end
c11=(1-fu11).*a1;
if trend==1
 c33=fu22.*cumsum(a1);
else
c33=fu22.*(a1);
end
c44=(1-fu22).*a1;

c1(isnan(pp))=0;c1(isnan(pp2))=0;c1(isnan(pp1))=0;c1(isnan(pp3))=0;
c2(isnan(pp))=0;c2(isnan(pp2))=0;c2(isnan(pp1))=0;c2(isnan(pp3))=0;
c11(isnan(pp))=0;c11(isnan(pp2))=0;c11(isnan(pp1))=0;c11(isnan(pp3))=0;
c22(isnan(pp))=0;c22(isnan(pp2))=0;c22(isnan(pp1))=0;a(isnan(pp))=0;a(isnan(pp2))=0;a(isnan(pp1))=0;a(isnan(pp3))=0;

if joint==1
    c1(isnan(pp3))=0;c2(isnan(pp3))=0;c11(isnan(pp3))=0;c22(isnan(pp3))=0;
end
% c1(a==0)=0;c2(a==0)=0;
c1=reshape(c1,size(c1,1),1,size(c1,2));c2=reshape(c2,size(c2,1),1,size(c2,2));%building cons for nonlinear
c11=reshape(c11,size(c11,1),1,size(c11,2));
if time_fe==0
c22=reshape(c22,size(c22,1),1,size(c22,2));%building cons for nonlinear
else
 c22=reshape(c22,size(c22,1),1,size(c22,2),size(c22,3));%building cons for nonlinear
end
c33=reshape(c33,size(c33,1),1,size(c33,2));c44=reshape(c44,size(c44,1),1,size(c44,2));%building cons for nonlinear
pp11=reshape(pp1,size(pp1,1),1,size(pp1,2));
pp55=reshape(pp5,size(pp5,1),1,size(pp5,2));


%   gzz=reshape(data2(2:end,12,:),size(data2,1)-1,size(data2,3));
if lev==0
gzz=t1(2:end,:);
else
gzz=t1(2:end,:);
end    
  %   gzz(1:end,:)=zscore(gzz(1:end,:));
  gzz(1:end,:)=1*(gzz(1:end,:));pp4=gzz(1:end-nlag,:);
 
clear gz
 for i=1:nlag
gz(:,:,i)=gzz(nlag+1-i:end-i+1,:);%gz variable
 end
%  size(tt),size(gz)
for i=1:nlag
     f=gz(:,:,i); 
     f(isnan(pp))=0; f(isnan(pp1))=0; f(isnan(pp2))=0;
 f(isnan(pp3))=0;
 if joint==1
     f(isnan(pp3))=0;
 end
 gz(:,:,i)=f;
 end
 for i=1:nlag-1
     f=y1(:,:,i);
 f(isnan(pp))=0;f(isnan(pp1))=0;f(isnan(pp2))=0;f(isnan(pp3))=0;
 if joint==1
     f(isnan(pp3))=0;
 end
 y1(:,:,i)=f;
 end
 y1(isnan(y1))=0;

for i=1:nlag-1
ful(:,:,i)=fu(nlag-i:end-i,:);%lagged dependent variable in first differences
end

%LINEAR
if time_fe==0
if near_var==0
parfor j=1:length(w)-1
  
  for i=1:hor
         if lev==0
     yy=tt(nlag+i:end,w(j):w(j+1)-1)- tt(nlag:end-i,w(j):w(j+1)-1); %variable entering in the LHS
         else
      yy=tt(nlag+i:end,:);
         end
           
gg=gz(1:end-i+1,w(j):w(j+1)-1,:);
 gg(gz(i:end,w(j):w(j+1)-1,:)==0)=0;
 if nlag>1
     if controls==1
gg1=y1(1:end-i+1,[w(j):w(j+1)-1,size(data3,2)+w(j):size(data3,2)+w(j+1)-1],:);
gg1(y1(i:end,[w(j):w(j+1)-1,size(data3,2)+w(j):size(data3,2)+w(j+1)-1],:)==0)=0;
     end
     if one_step==1
% gg1=y1(1:end-i+1,[w(j):w(j+1)-1,size(data3,2)+w(j):size(data3,2)+w(j+1)-1,...
%     2*size(data3,2)+w(j):2*size(data3,2)+w(j+1)-1,3*size(data3,2)+w(j):3*size(data3,2)+w(j+1)-1,...
%     4*size(data3,2)+w(j):4*size(data3,2)+w(j+1)-1,5*size(data3,2)+w(j):5*size(data3,2)+w(j+1)-1,...
%     6*size(data3,2)+w(j):size(data3,2)+w(j+1)-1,6*size(data3,2)+w(j):6*size(data3,2)+w(j+1)-1,...
%     7*size(data3,2)+w(j):7*size(data3,2)+w(j+1)-1],:);
% gg1(y1(i:end,[w(j):w(j+1)-1,size(data3,2)+w(j):size(data3,2)+w(j+1)-1,...
%     2*size(data3,2)+w(j):2*size(data3,2)+w(j+1)-1,3*size(data3,2)+w(j):3*size(data3,2)+w(j+1)-1,...
%     4*size(data3,2)+w(j):4*size(data3,2)+w(j+1)-1,5*size(data3,2)+w(j):5*size(data3,2)+w(j+1)-1,...
%     6*size(data3,2)+w(j):size(data3,2)+w(j+1)-1,6*size(data3,2)+w(j):6*size(data3,2)+w(j+1)-1,...
%     7*size(data3,2)+w(j):7*size(data3,2)+w(j+1)-1],:)==0)=0;
gg1=repelem(y1(1:end-i+1,1:8,:),1,w(j+1)-w(j));
gg1(repelem(y1(1:end-i+1,1:8,:),1,w(j+1)-w(j))==0)=0;

     end    
     
 else
gg1=y1(1:end-i+1,w(j):w(j+1)-1,:);
gg1(y1(i:end,w(j):w(j+1)-1,:)==0)=0;

 end
%  gg2222=fu(1:end-i+1,w(j):w(j+1)-1)+1-fu(1:end-i+1,w(j):w(j+1)-1);
%  gg2222((y(i:end,w(j):w(j+1)-1))==0)=0;
if nlag>1
    yy(gg(:,:,end)==0)=0;    
    if controls==1
      yy(gg1(:,1:w(j+1)-w(j),end)==0)=0;
            yy(gg1(:,w(j+1)-w(j)+1:2*(w(j+1)-w(j)),end)==0)=0;%one_Step you don't need to do this becasue of no missing obs for aggregate variables

    end  
%     if one_step==1
%       yy(gg1(:,1:w(j+1)-w(j),end)==0)=0;
%              yy(gg1(:,w(j+1)-w(j)+1:2*(w(j+1)-w(j)),end)==0)=0;%one_Step you don't need to do this becasue of no missing obs for aggregate variables
%              yy(gg1(:,2*(w(j+1)-w(j))+1:3*(w(j+1)-w(j)),end)==0)=0;%one_Step you don't need to do this becasue of no missing obs for aggregate variables
%              yy(gg1(:,3*(w(j+1)-w(j))+1:4*(w(j+1)-w(j)),end)==0)=0;%one_Step you don't need to do this becasue of no missing obs for aggregate variables
%              yy(gg1(:,4*(w(j+1)-w(j))+1:5*(w(j+1)-w(j)),end)==0)=0;%one_Step you don't need to do this becasue of no missing obs for aggregate variables
%              yy(gg1(:,5*(w(j+1)-w(j))+1:6*(w(j+1)-w(j)),end)==0)=0;%one_Step you don't need to do this becasue of no missing obs for aggregate variables
%              yy(gg1(:,6*(w(j+1)-w(j))+1:7*(w(j+1)-w(j)),end)==0)=0;%one_Step you don't need to do this becasue of no missing obs for aggregate variables
%              yy(gg1(:,7*(w(j+1)-w(j))+1:8*(w(j+1)-w(j)),end)==0)=0;%one_Step you don't need to do this becasue of no missing obs for aggregate variables
% 
%     end  
else
% yy(gg(:,1:-w(j)+w(j+1),end)==0)=0;
yy(gg(:,:,end)==0)=0;
end


aa=a(1:end-i+1,:,w(j):w(j+1)-1);  
aa(yy(1:end,:)==0)=0;
ct=num2cell(aa,[1 2]);

cons = blkdiag(ct{:});% basic code is a=rand(3,3,4);c=num2cell(a,[1 2]);M = blkdiag(c{:})




if gls==0
if treat==1
% x=[reshape(gg2222,numel(gg2222),1),...
%     cons{i},reshape(gg1,numel(gg1)/(nlag-1),nlag-1),...
%    reshape(gg,numel(gg)/nlag,nlag)];
else
    if nlag>1
if controls==1
    
    x1=cons.*repmat(reshape(gg(:,:,1),numel(gg(:,:,1))/1,1),1,size(cons,2));
    x1=sparse(x1);x2=cell(nlag-1,1);
    for z=1:nlag-1
%         x2=repmat(cons{i},1,2*nlag-2).*repmat(reshape(sparse(gg1(:,:,:)),size(cons{i},1),(2*nlag-2)),1,size(cons{i},2));
%      x2{1,z}=repmat(cons,1,2).*repmat(reshape(sparse(gg1(:,:,z)),size(cons,1),2),1,size(cons,2));%2 is number of controls
      x2{1,z}=repelem(cons,1,2).*repmat(reshape(sparse(gg1(:,:,z)),size(cons,1),2),1,size(cons,2));%2 is number of controls
%      save(['INPUT_V_X_XPXI_CONS/X2_BASELINE' num2str(z) '.mat'],'x2');

    end
    x2=cell2mat(x2);
      x2=sparse(x2);
   
x=[x1 x2 cons];
x=sparse(x);
end
if one_step==1

    x1=cons.*repmat(reshape(gg(:,:,1),numel(gg(:,:,1))/1,1),1,size(cons,2));
    x1=sparse(x1);x2=cell(nlag-1,1);
    for z=1:nlag-1
     x2{1,z}=repelem(cons,1,8).*repmat(reshape(sparse(gg1(:,:,z)),size(cons,1),8),1,size(cons,2));%8 is number of controls
%           x2{1,z}=repmat(cons,1,8).*repmat(reshape(sparse(gg1(:,:,z)),size(cons,1),8),1,size(cons,2));%2 is number of controls

    end

    x2=cell2mat(x2);
      x2=sparse(x2);
   
x=[x1 x2 cons];
x=sparse(x);%save('x','x');
end
    else
%    save('dd','cons','gg','yy'); 
if siz==0
if nonlin==0
%     if time_fe==0
       
            if fe==0%that means there are firm fe
x=[cons.*repmat(reshape(gg,numel(gg)/nlag,nlag),1,size(cons,2)),cons];
x=sparse(x);
            else
               
x=cons.*repmat(reshape(gg,numel(gg)/nlag,nlag),1,size(cons,2));
% x=sparse(x);
            end
    
else
 x=[cons.*repmat(reshape(gg,numel(gg)/nlag,nlag),1,size(cons,2)),...
    cons.*repmat(reshape(gg.^2,numel(gg)/nlag,nlag),1,size(cons,2)),...
    cons];
x=sparse(x); 
end
else
    x=[cons,reshape(gg.^3,numel(gg)/nlag,nlag),reshape(gg,numel(gg)/nlag,nlag),...
      reshape(gg.^2,numel(gg)/nlag,nlag)];
end
    end
end
else
    x=[cons{i},repmat(cons{i},1,nlag-1).*repmat(reshape(gg1,numel(gg1)/(nlag-1),nlag-1),1,size(data2,3)),...
   reshape(repmat(cons{i}',1,nlag)',size(cons{i},1),size(data2,3)*nlag).*repmat(reshape(gg,numel(gg)/(nlag),nlag),1,size(data2,3))];  
end

 yyy=reshape(yy,numel(yy),1);
 %yy(isnan(yy))=0;yy(yy==Inf)=0;
% save('s','yyy','x','pp','pp1','pp2','pp3','pp4');
    results=nwest_balanced_panel(yyy,x,i,size(yy,1),size(yy,2),nlag,no_interaction,gls,rep,time_fe,fe);
    clear_x;  
    b=results.beta;
%     bint=[b-(results.se*sig) b-(results.se*sig)];
%     bint(:,2)=b+(results.se*sig);

    reg{j}(:,i)=b(1:w(j+1)-w(j),:); 
    if nonlin==1
       reg1{j}(:,i)=b(w(j+1)-w(j)+1:(w(j+1)-w(j))*2,:); 
    end
 end
 end
if nonlin==0
reg=vertcat(reg{:});
else
reg=vertcat(reg{:});
reg1=vertcat(reg1{:});
reg=[reg;reg1];

end
else
  load data_baseline_with_rpi;
data(:,1)=[data(1,1);diff(data(:,1))];
data(:,3)=[];data(:,end)=[];


    cons=zeros(1,size(data2,3));
rr=data2(2:end,or,:);
rr(isnan(rr))=0;
parfor i=1:size(data2,3)
    s=rr(:,:,i);
     [~,~,usa]=bootstrap_ISTC_news([data(:,:) s],1,5,20,20,1,3,1,1,1,1,size(data,2)-1,0);
  reg(i,:)=usa.TFP(:,end)'; 
end

end
else
    %START TIME_FE ESTIMATION
       u=eye((size(data2,1)-1));
 parfor j=1:length(w)-1
  for i=1:hor
         if lev==0
     yy=tt(nlag+i:end,w(j):w(j+1)-1)- tt(nlag:end-i,w(j):w(j+1)-1); %variable entering in the LHS
         else
      yy=tt(nlag+i:end,:);
         end
           
gg=gz(1:end-i+1,w(j):w(j+1)-1,:);
gg(gz(i:end,w(j):w(j+1)-1,:)==0)=0;
gg1=y1(1:end-i+1,w(j):w(j+1)-1,:);
gg1(y1(i:end,w(j):w(j+1)-1)==0,:)=0;
%  gg2222=fu(1:end-i+1,w(j):w(j+1)-1)+1-fu(1:end-i+1,w(j):w(j+1)-1);
%  gg2222((y(i:end,w(j):w(j+1)-1))==0)=0;

if nlag>1
      yy(gg1(:,w(j):w(j+1)-1,end)==0)=0;
else
yy(gg(:,:,end)==0)=0;
end

aa=a(1:end-i+1,:,w(j):w(j+1)-1);  
aa(yy(1:end,1:-w(j)+w(j+1))==0)=0;
ct=num2cell(aa,[1 2]);

cons = blkdiag(ct{:});% basic code is a=rand(3,3,4);c=num2cell(a,[1 2]);M = blkdiag(c{:})

% for q=1:size(u,2)
%    uu=repmat(u(1:end-i+1,q),w(j+1)-w(j),1);
%  uu(yy(1:end,1:-w(j)+w(j+1))==0)=0;
% uu1{i,j}(:,q)=circshift([uu;zeros(size(cons{i},1)-size(uu,1),1)],(w(j+1)-w(j))*size(uu,1));
 uu=repmat(u(1:end-i+1,:),w(j+1)-w(j),1);
 z=reshape(yy(:,1:-w(j)+w(j+1)),(-w(j)+w(j+1))*size(yy,1),1);
 z=repmat(z,1,size(uu,2));
uu(z==0)=0;
% end
%   uu1{i,j}=reshape(uu1{i,j},size(uu1{i,j},1),size(uu1{i,j},2)*size(uu1{i,j},3));     
  x=[cons.*repmat(reshape(gg,numel(gg)/nlag,nlag),1,size(cons,2)),...
      cons,uu];
  %x(isnan(x))=0;x(x==Inf)=0;
  
% x=sparse(x);

 yyy=reshape(yy,numel(yy),1);
 %yy(isnan(yy))=0;yy(yy==Inf)=0;
% save('s','yyy','x','pp','pp1','pp2','pp3','pp4');
    results=nwest_balanced_panel(yyy,x,i,size(yy,1),size(yy,2),nlag,no_interaction,gls,rep,time_fe,fe);
    clear_x;  
    b=results.beta;
%     bint=[b-(results.se*sig) b-(results.se*sig)];
%     bint(:,2)=b+(results.se*sig);

    reg{j}(:,i)=b(1:w(j+1)-w(j),:);
    
 end
 end
reg=vertcat(reg{:});

end
   reg_firms=reg;        

 if nonlin==0
    cons=zeros(1,size(reg,1));

reg_lin=median(reg_firms);%MEDIAN MPK FIRM RESPONSE
% reg_low=median(reg(1:round(quantile(1:size(cons,2),0.25)),:));%MEDIAN MPK FIRM RESPONSE
% reg_up=median(reg(round(quantile(1:size(cons,2),0.75)):size(cons,2),:));%MEDIAN MPK FIRM RESPONSE
if near_var==0
mis=data3*reg_firms;
else
 index1 = find(mpk>quantile(mpk,0.30), 1, 'first');
index2 = find(mpk>quantile(mpk,0.70), 1, 'first'); 
mis=data3(:,index1:index2)*reg_firms(index1:index2,:,:);
end
index1 = find(mpk>quantile(mpk,0.25), 1, 'first');
index2 = find(mpk>quantile(mpk,0.75), 1, 'first');

% index3 = find(size1>quantile(size1,0.25), 1, 'first');
% index4 = find(size1>quantile(size1,0.75), 1, 'first');
% index5 = find(cap_share>quantile(cap_share,0.25), 1, 'first');
% index6 = find(cap_share>quantile(cap_share,0.75), 1, 'first');
% [data3,all_ord]=sort(data3);
% 
% index7 = find(data3>quantile(data3,0.25), 1, 'first');
% index8 = find(data3>quantile(data3,0.75), 1, 'first');
% if time_fe==1
reg(1:size(cons,2),:)=reg(mpk_ord,:);
% end
reg_low=median(reg(1:index1,:));%MEDIAN MPK FIRM RESPONSE
reg_up=median(reg(index2:size(cons,2),:));%MEDIAN MPK FIRM RESPONSE
reg(1:size(cons,2),:)=reg(size_ord,:);
reg_low_size=median(reg(1:index1,:));%MEDIAN MPK FIRM RESPONSE
reg_up_size=median(reg(index2:size(cons,2),:));%MEDIAN MPK FIRM RESPONSE
reg(1:size(cons,2),:)=reg(cap_share_ord,:);
reg_low_cap=median(reg(1:index1,:));%MEDIAN MPK FIRM RESPONSE
reg_up_cap=median(reg(index2:size(cons,2),:));%MEDIAN MPK FIRM RESPONSE
reg(1:size(cons,2),:)=reg(all_ord,:);
reg_low_all=median(reg(1:index1,:));%MEDIAN MPK FIRM RESPONSE
reg_up_all=median(reg(index2:size(cons,2),:));%MEDIAN MPK FIRM RESPONSE

reg_asym=reg_up-reg_low;%mpk response difference
reg_asym_size=reg_up_size-reg_low_size;%size response difference
reg_asym_cap=reg_up_cap-reg_low_cap;%size response difference
reg_asym_all=reg_up_all-reg_low_all;%size response difference

reg_lin1=reg_lin;mis1=mis;reg_low1=reg_low;reg_up1=reg_up;
reg_asym1=reg_asym;reg_low_size1=reg_low_size;
reg_up_size1=reg_up_size;reg_asym_size1=reg_asym_size;
    reg_low_cap1=reg_low_cap;reg_up_cap1=reg_up_cap;
    reg_asym_cap1=reg_asym_cap;reg_low_all1=reg_low_all;
    reg_up_all1=reg_up_all;reg_asym_all1=reg_asym_all;
 reg_lin2=reg_lin;mis2=mis;reg_low2=reg_low;reg_up2=reg_up;
reg_asym2=reg_asym;reg_low_size2=reg_low_size;
reg_up_size2=reg_up_size;reg_asym_size2=reg_asym_size;
    reg_low_cap2=reg_low_cap;reg_up_cap2=reg_up_cap;
    reg_asym_cap2=reg_asym_cap;reg_low_all2=reg_low_all;
    reg_up_all2=reg_up_all;reg_asym_all2=reg_asym_all;
      
    
 else
         cons=zeros(1,size(data3,2));

     %POSITIVE SHOCK IRF
reg_firms=reg(1:size(cons,2),:)+reg(size(cons,2)+1:2*size(cons,2),:);
reg_lin=median(reg_firms);%MEDIAN MPK FIRM RESPONSE
mis=data3*reg_firms;

index1 = find(mpk>quantile(mpk,0.25), 1, 'first');
index2 = find(mpk>quantile(mpk,0.75), 1, 'first');


reg_low=median(reg_firms(1:index1,:));%MEDIAN MPK FIRM RESPONSE
reg_up=median(reg_firms(index2:size(cons,2),:));%MEDIAN MPK FIRM RESPONSE
reg_firms(1:size(cons,2),:)=reg_firms(size_ord,:);
reg_low_size=median(reg_firms(1:index1,:));%MEDIAN MPK FIRM RESPONSE
reg_up_size=median(reg_firms(index2:size(cons,2),:));%MEDIAN MPK FIRM RESPONSE
reg_firms(1:size(cons,2),:)=reg_firms(cap_share_ord,:);
reg_low_cap=median(reg_firms(1:index1,:));%MEDIAN MPK FIRM RESPONSE
reg_up_cap=median(reg_firms(index2:size(cons,2),:));%MEDIAN MPK FIRM RESPONSE
reg_firms(1:size(cons,2),:)=reg_firms(all_ord,:);
reg_low_all=median(reg_firms(1:index1,:));%MEDIAN MPK FIRM RESPONSE
reg_up_all=median(reg_firms(index2:size(cons,2),:));%MEDIAN MPK FIRM RESPONSE

reg_asym=reg_up-reg_low;%mpk response difference
reg_asym_size=reg_up_size-reg_low_size;%size response difference
reg_asym_cap=reg_up_cap-reg_low_cap;%size response difference
reg_asym_all=reg_up_all-reg_low_all;%size response difference
     
     %NEGATIVE SHOCK IRF
     reg_firms1=reg_firms;
reg_firms=reg(1:size(cons,2),:)-reg(size(cons,2)+1:2*size(cons,2),:);
reg_lin1=median(reg_firms);%MEDIAN MPK FIRM RESPONSE
mis1=data3*reg_firms;

index1 = find(mpk>quantile(mpk,0.25), 1, 'first');
index2 = find(mpk>quantile(mpk,0.75), 1, 'first');


reg_low1=median(reg_firms(1:index1,:));%MEDIAN MPK FIRM RESPONSE
reg_up1=median(reg_firms(index2:size(cons,2),:));%MEDIAN MPK FIRM RESPONSE
reg_firms(1:size(cons,2),:)=reg_firms(size_ord,:);
reg_low_size1=median(reg_firms(1:index1,:));%MEDIAN MPK FIRM RESPONSE
reg_up_size1=median(reg_firms(index2:size(cons,2),:));%MEDIAN MPK FIRM RESPONSE
reg_firms(1:size(cons,2),:)=reg_firms(cap_share_ord,:);
reg_low_cap1=median(reg_firms(1:index1,:));%MEDIAN MPK FIRM RESPONSE
reg_up_cap1=median(reg_firms(index2:size(cons,2),:));%MEDIAN MPK FIRM RESPONSE
reg_firms(1:size(cons,2),:)=reg_firms(all_ord,:);
reg_low_all1=median(reg_firms(1:index1,:));%MEDIAN MPK FIRM RESPONSE
reg_up_all1=median(reg_firms(index2:size(cons,2),:));%MEDIAN MPK FIRM RESPONSE

reg_asym1=reg_up1-reg_low1;%mpk response difference
reg_asym_size1=reg_up_size1-reg_low_size1;%size response difference
reg_asym_cap1=reg_up_cap1-reg_low_cap1;%size response difference
reg_asym_all1=reg_up_all1-reg_low_all1;%size response difference
       
     %POSTIIVE MINUS NEGATIVE SHOCK IRF
     reg_firms2=reg_firms;
reg_firms=2*reg(size(cons,2)+1:2*size(cons,2),:);
reg_lin2=median(reg_firms);%MEDIAN MPK FIRM RESPONSE
mis2=data3*reg_firms;

index1 = find(mpk>quantile(mpk,0.25), 1, 'first');
index2 = find(mpk>quantile(mpk,0.75), 1, 'first');


reg_low2=median(reg_firms(1:index1,:));%MEDIAN MPK FIRM RESPONSE
reg_up2=median(reg_firms(index2:size(cons,2),:));%MEDIAN MPK FIRM RESPONSE
reg_firms(1:size(cons,2),:)=reg_firms(size_ord,:);
reg_low_size2=median(reg_firms(1:index1,:));%MEDIAN MPK FIRM RESPONSE
reg_up_size2=median(reg_firms(index2:size(cons,2),:));%MEDIAN MPK FIRM RESPONSE
reg_firms(1:size(cons,2),:)=reg_firms(cap_share_ord,:);
reg_low_cap2=median(reg_firms(1:index1,:));%MEDIAN MPK FIRM RESPONSE
reg_up_cap2=median(reg_firms(index2:size(cons,2),:));%MEDIAN MPK FIRM RESPONSE
reg_firms(1:size(cons,2),:)=reg_firms(all_ord,:);
reg_low_all2=median(reg_firms(1:index1,:));%MEDIAN MPK FIRM RESPONSE
reg_up_all2=median(reg_firms(index2:size(cons,2),:));%MEDIAN MPK FIRM RESPONSE

reg_asym2=reg_up2-reg_low2;%mpk response difference
reg_asym_size2=reg_up_size2-reg_low_size2;%size response difference
reg_asym_cap2=reg_up_cap2-reg_low_cap2;%size response difference
reg_asym_all2=reg_up_all2-reg_low_all2;%size response difference
 reg_firms3=reg_firms;
 reg_firms=[reg_firms1;reg_firms2;reg_firms3];%POSITIVE;NEGATIVE;ASYMMETRY
     
 end
 
clear c* g* a* f*
